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ABSTRACT 

We apply standard Markov chain Monte Carlo (MCMC) analysis techniques for 50 short- 
period, single-planet systems discovered with radial velocity technique. We develop a new 
method for accessing the significance of a non-zero orbital eccentricity, namely F analysis, 
which combines frequentist bootstrap approach with Bayesian analysis of each simulated data 
set. We find the eccentricity estimations from F analysis are generally consistent with results 
from both standard MCMC analysis and previous references. The T method is particular use- 
ful for assessing the significance of small eccentricities. Our results suggest that the current 
sample size is insufficient to draw robust conclusions about the roles of tidal interaction and 
perturbations in shaping the eccentricity distribution of short-period single-planet systems. 
We use a Bayesian population analysis to show that a mixture of analytical distributions is a 
good approximation of the underlying eccentricity distribution. For short-period planets, we 
find the most probable values of parameters in the analytical functions given the observed 
eccentricities. These analytical functions can be used in theoretical investigations or as pri- 
ors for the eccentricity distribution when analyzing short-period planets. As the measurement 
precision improves and sample size increases, the method can be applied to more complex 
parametrizations for the underlying distribution of eccentricity for extrasolar planetary sys- 
tems. 
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1 INTRODUCTION 

The discovery of exoplanets has significantly advanced our 
unde rstanding of formation and evolution o f planetary sys- 
tem jWolszczanl 1 1994 iMavor & Oueloj 1 19951 : iMarcv & Butler] 
1 19961) . As of February 20 1 1 , over 500 exoplanets have been discov- 
ered including 410 systems detected by radial velocity (RV) tech- 
nique 1 . The eccentricity distribution of exoplanets is very different 
from that of solar system. For sufficiently short-period planets, it 
is expected that tidal circularization would lead to nearly circular 
orbits. Yet, several short-period planets appear to have eccentric 
orbits. Several mechanisms (e.g. planet scattering, Kozai effect) 
have been proposed to explain the observed eccentricity distribu- 
tion jTakeda & Rasid20o"5r . lZhou & Linl2007l : lFord & Rasiol2008l : 



ljuric & Tremaind 20081) . This paper aims to improve our under- 
standing of the true eccentricity distribution and its implications 
for orbital evolution. 

The Bayesian approach offers a rigorous basis for determin- 
ing the posterior eccentricity distribution for individual system. 
The Bayesian method is particularly advantageous relative to tra- 
ditional bootstrap method when the orbital eccentricity is poorly 
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constrained by RV data dFordl2006l) . [Fot3 J2006t) discussed eccen- 
tricity estimation using Markov Chain Monte Carlo (MCMC) sim- 
ulation in the framework of Bayesian inference theory and found a 
parameter set that accelerates convergence of MCMC for low ec- 
centricity orbit. For a population of planets on nearly circular orbits, 
eccentricity estimates for planets on circular orb it are biased re- 
sultin g in overestimation of orbital eccentricities dZakamska et al.l 
l201ll) . Further complicating matters, the population of known exo- 
planets is not homogeneous, and the observed eccentricity distribu- 
tion is affected by the discovery method, selection effects and data 
analysis technique. 

In this paper, we construct a catalog of short-period single- 
planet systems using homogeneous RV data reduction process, 
i.e., standard MCMC analysis (§2.1). In §2.2, we describe a new 
method of estimating Keplerian orbital parameters, namely T anal- 
ysis. We present the results for standard MCMC analysis in §3. 
We compare the results of both methods with each other and the 
results from previous references. We interpret the results to investi- 
gate how tidal effect and perturbation affect the orbital eccentricity 
distribution of short-period planets (§4). In §5, we provide an an- 
alytical function for the underlying eccentricity distribution that is 
able to reproduce the observed the eccentricity distribution. We dis- 
cuss the results in §6 and summarize our conclusions in §7. 
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2 METHOD 

We select all the systems with: 1) a single known planet discovered 
with the radial velocity technique as of April 2010; 2) an orbital pe- 
riod of less than 50 days; and 3) a publicly available radial velocity 
data set. We exclude planets discovered by the transit technique in 
order to avoid complications due to selection effects dGaudi et al.l 
l2005h . We perform an orbital analysis on each system in our sample 
using: 1) a standard MCMC analysis (§2.1) and 2) a new method, 
T analysis (described in §2.2). We focus on the eccentricity estima- 
tion for each planet since the eccentricity is an important indication 
of orbital evolution and tidal interaction. 



2.1 Bayesian Orbital Analysis of Individual Planet 

We performed a Bayesian analysis of the published radial velocity 
observation using a model consisting of one low-mass companion 
following a Keplerian orbit. If a long-term RV trend is included 
in the original paper reporting the RV data or if a linear trend of 
more than 1 m • s • yr~' is apparent, then we add to the model a 
constant long-term acceleration due to distant planetary or stellar 
companion. 

We calculate a posterior sample using the Markov Chai n 
Monte Carlo (MCMC) technique as described in iFord d2006l) . 
Each state in the Markov chain is described by the parameter set 
6 = {P,K,e,aj,M ,C h D,<Tj}, where P is orbital period, K is the 
velocity semi-amplitude, e is the orbital eccentricity, a> is the argu- 
ment of periastron, M is the mean anomaly at chosen epoch r, C, is 
constant velocity offset (subscript i indicates constant for different 
observatory), D is the slope of a long-term linear velocity trend, and 
CTj is the "jitter" parameter. The jitter parameter describes any addi- 
tional noise includin g both astrophysica l noises, e.g., stellar oscil- 
lation, stellar spots dWright et al .1120050 and any instrument noise 
not accounted for in the reported measurement uncertainties. The 
RV perturbation of a host star at time t k due to a planet on Keplerian 
orbit and possible perturbation is given by 

v k g = K ■ [cos(w + T) + e ■ cosfw)] + D ■ (t k - t), (1) 

where T is the true anomaly which is related to eccentric anomaly 
E via the relation, 

tan 9= Vr^H!)- (2) 

The eccentric anomaly is related to the mean anomaly M via Ke- 
pler's equation, 

E(t) - e ■ sin[E(f)] = M(t) - M = y{t - t). (3) 

We choose prio rs of each parameter as described in 
lFord& Gregory] d2007t) . The prior is uniform in logarithm of orbital 
period. For K and < x, we use a mo dified Jefferys prior in the form of 
p(x) oc (x + Xo)' 1 jGregorvil2005h with K mi „ = o~j min = 0.1 m ■ s" 1 . 
Priors are uniform for: e (0 < e < 1), a> and M (0 < u>, M < 2n), 
C, and D. We verified that the parameters in 6 did not approach the 
limiting values. We assume each observation results in a measure- 
ment drawn from normal distribution centered at the true velocity, 
resulting in a likelihood (i.e., conditional probability of making the 
specified measurements given a particular set of model parameters) 
of 

p-r exp[-(y kS - v k ) 2 /2a-j] 
pW, M) oc [I £» — , (4) 





80 












HD_68988 : 




60 
















40 










68% 




Dens! 


20 











* 

o 


v . X ■ 


t_ 

<u 
n 


80 




0.05 


0.10 


0.15 


E 

3 












HD_330075 


Z 


60 

40 

20 



68% 

* 

o 


<t) 
Q. 
Q. 
D 











0.1 0.2 0.3 0.4 

e 



Figure 1. Examples of how credible intervals of standard MCMC analysis 
are calculated using posterior distribution of e. Grey region contains 68% 
of total number of posterior samples of e. 



where v k is radial velocity at time t k , and v ke is the model velocity 
at time t k given the model parameters 6. Noise cr k consists of two 
parts. One component is from the observation uncertainty a k ^ s re- 
ported in the radial velocity data, and the other is the jitter, tr ; , 
which accounts for any unforseen additional noise including in- 
strument instability and stellar jitter. The two parts are added in 
quadrature in order to generate a k . We calculate the Gelman-Rubin 
statistic, R, to test for nonconvergence of Markov chains. 

We perform a MCMC analysis for RV data set of each sys- 
tem in the sample and obtain posterior samples of h and k, where 
h = e cos a> and k = e sin oj. This parameterization has been shown 
to be more effective in description of the eccentricity distribution 
for low eccentric orbits dFord]|2006l) . We take steps in h and k and 
adjust the acceptance rate according to the Jacobian of the coor- 
diante transformation, so as to maintain a prior that is uniform in e 
and at. Mean values, h and k, from posterior samples of h and k are 
adopted to calculate e^cmc using the equation €mcmc = 

vPTF. 

The posterior distribution of e is not always Gaussian distribution 
especially near e ~0. Therefore, it is not appropriate to calculate the 
uncertainty of e using the equation of error propagation in which 
gaussian noise is assumed. We use posterior distribution of e to in- 
fer the credible interval of e. The boundaries of the region where 
68% posterior samples populate are adopted as e\ mter and e upper 
(Appendix I A It . Fig. Q] illustrates two examples of how the credi- 
ble intervals are inferred for HD 68988 (eccentric orbit) and HD 
330075 (circular orbit). 



2.2 T Analysis of Individual Systems 

A fully Bayesian analysis of the population of exoplanet eccen- 
tricities would be computationally prohibitive due to the large 
number of dimensions. Therefore, we develop a hybrid Bayesian- 
frequentist method to evaluate the significance of a non-zero ec- 
centricity measurement. We combine a bootstrap style approach of 
generating and analyzing synthetic data sets with MCMC analysis 
of each synthetic data set to obtain a frequentist confidence level 
for each eccentricity that accounts for biases. First, we perform the 
standard MCMC analysis described in §2.1 on the real RV data 
set and adopt the mean value of each orbital parameter in 8 ex- 
cept e. We generate a series of simulated radial velocity data sets 
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Figure 2. Contours of posterior distribution in h and k space for HD 68988 
(Solid-68% of sample points included; dashed-95% of sample points in- 
cluded; dotted-99% of sample points included). Red contours are poste- 
rior distribution for real observation, green contours are for simulated RV 
data set with e=0.00, and blue contours are for simulated RV data set 
with e=0.13 Ecce ntricity of HD 68988 is 0.1250±0.0087 according to 
iButler et alj 120061) . 



at different values of e. The adopted K is scaled accordingly to 
K oc (1 - e 2 y - 5 . The simulated radial velocity data has the same 
number of observations, and each simulated observation takes place 
at exactly the same time and the same mean anomaly as the real ob- 
servation. Gaussian noise with standard deviation of at (§2.1) are 
added to simulated radial velocity data sets at different eccentrici- 
ties. Each simulated RV data set has the same reported RV measure- 
ment uncertainties as the real RV observations. Standard MCMC 
analysis is then performed on each of the simulated RV data sets. 

For both real and simulated data sets, we construct a two- 
dimensional histogram using the posterior samples in (h, k) space 
to approximate a two-dimensional posterior distribution for h 
and k, d,{hj,kj) and d s (hj,kj), where i and j denote bin in- 
dices, and the subscripts r and s denote the real and sim- 
ulated data set. We compare the distribution for each simu- 
lated data set d s (hj,kj) to the distribution for real radial ve- 
locity data set d,.(hj,kj). To quantify the similarity between 
d s (h,,kj) and d r (hi,kj), we calculate the statistic defined as T = 

N N 

iVV( 



^ jJAJQh, kj) - ft,) ■ (d r (h„ kj) - ^i r y\l[a s a r (N 2 - 1)], where N 

i=l 7=1 

is the number of bins in h or k dimension, \i and <x represent mean 
and standard deviation. In other words, the T statistic is obtained 
by cross-correlating two posterior distributions in h and k space. 
Fig [2] illustrates the process by which we obtain T for the case of 
HD 68988. If d s {h„kj) matches d r {hi,kj), we expect to obtain a T 
value that approaches unity (blue and red contours). If the samples 
differ significantly then T decreases towards zero (red and green 
contours). For each eccentricity, we simulated 21 radial velocity 
data sets and compare d s (h h kj) with d r (h h kj) to obtain 21 T statis- 
tics between simulated and real RV data. We choose the median 
value f as an indicator of overall similarity at given eccentricity. 

Based on above analysis of simulated radial velocity data with 
different input eccentricities e, we obtain a relationship between T 
and e, i.e. T(e). We use a high-order polynomial to interpolate for 
f (e). We define e to be the eccentricity at which f (e) reaches max- 
imum and we interpret e as an estimator of eccentricity. For HD 
68988 (Fig. [3), input eccentricity ranges from 0.00 to 0.29 with 
step size of 0.01. f (e) reaches maximum at e = 0.134. We estimate 



Figure 3. f as a function of eccentricity e for HD 68988. Open circles are 
results from simulations, solid line is the result of polynomial fitting. The 
long dashed line is the critical threshold, r c o.68 at the 68% confidence level. 



statistical confidence level of e using every pair of posterior eccen- 
tricity samples calculated from the data sets that are generated as- 
suming the same eccentricity. Consider the example of HD 68988 
again: 21 sets of posterior distribution in h and k space, d s (hj,kj), 
are obtained. Comparison between each pair gives a T statistic be- 
tween simulated RV data. 2)?:, = 210 T statistics in total for simu- 
lated RV data sets are calculated at eccentricity of 0.13 and 68.1% 
of pairs (143 out of 210) have a T statistic greater than 0.3714. We 
define this value, r C] o.68. as me critical T value for HD 68988 at 68% 
confidence level for e = 0.13. Therefore, if T statistic obtained in 
comparison between d s {hj,kj) and d r (hj,kj) is less than 1^0.68 > we 
argue that the eccentricity inferred from simulated RV data set is 
not consistent with the observed eccentricity of the system at 68% 
confidence level. In the case where e is located between grids of 
simulated e values, we calculate T Ct o.6s at e using interpolation of 
r t -,0 68 at nearby e values. We use a high-order polynomial to ap- 
proximate the discrete data T(e). The polynomial is later used to 
infer e, lower and upper limit of eccentricity. For HD 68988, r Ci o.68 
is 0.3663 at e = 0.134 after interpolation. Using the relationship 
between f and e (Fig. [3}, we look for the e values corresponding 
to r Cj o,68 as estimators of the lower and upper limit for eccentric- 
ity of the planet system at a 68% confidence level. In HD 68988, 
we obtained e = 0.1 34^ {J* using T analysis. In comparison, we 
have obtained e = 0. 119^n'mi using a standard MCMC analysis 



and IButler et alj J2006I) reported e = 0.125 + 0.009. 



3 RESULTS FOR INDIVIDUAL PLANETS 

In our sample of 50 short-period single-planet systems, we success- 
fully analyzed 42 systems using T analysis, and 46 systems using 
standard MCMC analysis (Appendix IAU . All the error are based 
on a 68% confidence level (T analysis) or a 68% credible inter- 
val (MCMC). The unsuccessful cases in standard MCMC analy- 
sis are HD 189733, HD 219828, HD 102195 and GJ 176. In HD 
189733, most of the RV data points were taken during observation 
of Rossiter-McLaughlin effect, which is not modeled here. MCMC 
analysis fails for HD 219828 and GJ 176 because of low signal to 
noise ratio (K=7 m ■ s~' and n obs =20 for HD 219828 and K=4.1 
m • s~ l and n obs =57 for GJ 176). RV data points of HD 102195 
were taken at 3 observatories and MCMC analysis is complicated 
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Figure 4. Top: comparison between standard MCMC analysis (blue) and 
previous references (red); Bottom: comparison between standard MCMC 
analysis (blue) and T analysis (green). The systems where disagreements 
take place are marked with a number: 1, HD 149026; 2, r Boo; 3, HD 
195019. 



using different observatory offsets. In addition to the above, T anal- 
ysis was unsuccessful for HD 162020, GJ 86, HD 17156 and HD 
6434. Since T analysis involves generating simulated RV data, the 
limited number of observations and partial phase coverage for these 
systems can cause poor convergence for some simulated data sets. 
These limitations become more severe for systems with high ec- 
centricity (e.g. HD 17156, e = 0.684) since phase coverage is more 
important for eccentric orbits. 



3.1 Comparison: Standard MCMC and References 

In Fig. [4] (top), we compare results from two sources, standard 
MCMC analysis and previous references. In most cases the two 
methods provide similar results. We found there are 3 systems for 
which the eccentricity estimates are not consistent, i.e., the pub- 
lished eccentricity error bar does not overlap the 68% credible in- 
terval from our MCMC an alysis. They are H P 149026, r Boo, HD 
195019. For HD 149026, ISato et al.1 J2005t) set the eccentricity to 
be zero when fitting the orbit. In contrast, we treat eccentricity as a 
variable and the standard MCMC method found that the 68% credi- 
ble interval for the systems mentioned above does not include zero. 
In addition, HD 149026 b is a known transiting planet and transit 
photometry provides additional constraint s on eccentricity which 
we have not included jCharbonnea3l2003l) . llCnutson et ail J2009h 
measured At n = 20.9^^ minute(2.9cr) for HD 149026 which is 
inconsistent with zero eccentricity, because e > ecosw ^ ^Af n , 
where P is period and Af n is the deviation of secondary eclipse 
from midpoint of primary transits. Standard MCMC results for 
other two systems (r Boo and HD 195019) are not consistent 
with those from previous references eve n though eccentrici ty was 
treated as variable in previous references. iButler et al.l |2006) report 
e = 0.023 ± 0.015 for t Boo and e = 0.014 ± 0.004 for HD 195019. 
On the contrary, standard MCMC analysis gives e=0.0787+[j ^ for 
t Boo and e=0.0017+S "^ for HD 195019 (See Appendix |7\T>. 



3.2 Comparison: Standard MCMC and T Analysis 

Fig. E] (bottom) compares the results from standard MCMC anal- 
ysis and T analysis. We find that the 68% credible/confidence in- 
tervals for the two methods overlap in the cases where there are 



discrepancies between standard MCMC analysis and previous ref- 
erences (see § 13. It . The confidence interval from the T analysis is 
generally larger than the credible interval from a standard MCMC 
analysis. The larger uncertainty from T analysis is likely due to the 
analysis accounting for the uncertainty in each velocity observation 
twice, first when generating synthetic data sets and a second time 
when analyzing the simulated data. Thus, the T analysis results in 
slightly larger uncertainty in eccentricity estimation. 

In order to understand the behavior of standard MCMC anal- 
ysis and T analysis for planets on nearly circular orbits, we con- 
duct an additional experiment generating many synthetic data sets 
where each system has a single planet on a circular orbit. We as- 
sume that they are observed at the same times and with the same RV 
measurement precisions as actual RV data sets. In order to under- 
stand the bias of each method for nearly circular systems, we com- 
pare the output eccentricities and their uncertainties. Using stan- 
dard MCMC analysis, we find that 76.4 ± 2.9% of the simulated 
data sets are consistent with zero using a 68% credible interval, 
and 23.6 ± 1.6% of the simulated data sets have 68% credible in- 
tervals that do not include zero. In contrast, for 16.8 ± 1.4% of the 
simulated data sets, the T analysis does not result in a 68% confi- 
dence interval that includes zero. In both cases, more than 68% of 
data sets are consistent with a circular orbit at a 68% level using ei- 
ther method. Using the T analysis, 6.8 ± 3.0% more simulated data 
sets are consistent with a circular orbit than based on the standard 
MCMC analysis. This confirms our intuition that the T analysis is 
a less biased method for analyzing systems at very small eccen- 
tricity. Thus, the T analysis may be a useful tool in assessing the 
significance of a measurement of a small non-zero eccentricity. In 
particular, we find 5 cases (1 1.4%) in which e; ower =0 for T analysis 
even though ei m . er for MCMC is greater than zero (e.g., HD 46375, 
HD 76700, HD 7924, HD 168746, HD 102117). 

A similar experiment is conducted except that an eccentric- 
ity of 0.2 is assigned to each system instead of zero eccentricity. 
Again, we assess the accuracy of the two methods by comparing 
the input and output eccentricities. When using the MCMC method, 
we find that the 68% credible interval for the eccentricity does not 
include the input eccentricity for 26.0 ± 2.1% of simulated data 
sets. When using the T method, we find that the 68% confidence 
interval for the eccentricity does not include the input eccentric- 
ity for 18.9 ± 1.8% of simulated data sets. Again, there is a larger 
fraction of results from the MCMC method that are not consistent 
with the input at a sizable eccentricity, indicating T analysis is less 
likely to reject the correct eccentricity. We also investigate the bias 
of the two methods at a significant eccentricity (i.e. e=0.2). In the 
cases where the output eccentricities are consistent with the input, 
we find that 47.9 ± 3.3% of the output eccentricities are below 0.2 
while 52.1 ± 3.5% of outputs are over 0.2 for MCMC method in- 
dicating the MCMC method is no t biased at a sizable ecc entricity, 
which agrees with the finding from lZakamska et al.l J201 lh . In com- 
parison, T analysis is also an unbiased analysis with 49.9 ± 3.3% 
below input and 50.1 ± 3.3% exceeding input. Therefore, we find 
no evidence for significant bias of either method for data sets with 
a significant eccentricity. 

3.3 Discussion of T Analysis 

The different methods for analyzing Doppler observations are com- 
plimentary. Bayesian methods and MCMC in particular are rou- 
tinely used to sample from the posterior distribution for the Ke- 
plerian orbital parameters for a given system. However, the anal- 
ysis of an exoplanet population is more complicated than simply 
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performing a Bayesian analysis of each system. To illustrate this 
point, consider a population of planets that all have exactly circu- 
lar orbits. Due to measurement uncertainties and finite sampling, 
the "best-fit" eccentricity for each system will be non-zero. Simi- 
larly, since eccentricity is a positive-definite quantity, the analysis 
of each system will result in a posterior distribution that has signifi- 
cant support for e > 0. This property remains even if one combines 
many point estimates (e.g., "best-fit"), frequentist confidence inter- 
vals or Bayesian posterior distributions. While the posterior distri- 
bution for the orbital parameters represents the best possible anal- 
ysis of an individual system, the inevitable bias for nearly circular 
orbits is a potential concern for population analyses. Therefore, it is 
important to apply different methods for p opulation analyses (e.g., 
iHogg et alj l l2O10l) ; IZakamska et al.l l l201lh ). 

Since we intend to investigate the potential role of tidal ef- 
fects on the eccentricity distribution of short-period planets, we 
developed a hybrid technique to assess the sensitivity of our re- 
sults to bias in the posterior distribution for planets with nearly cir- 
cular orbits. This hybrid technique (r analysis) involves perform- 
ing Bayesian analyses of each individual planetary system along 
with several simulated data sets, each of varying eccentricity. The 
MCMC analysis of each data set allows us to account for the vary- 
ing precision of eccentricity measurements depending on the ve- 
locity amplitude, measurement precision, number of observations 
and phase coverage. We assess the extent of the eccentricity bi- 
ases by performing the same analysis on simulated data sets with 
known eccentricity. We compare the posterior densities for the ac- 
tual data set to the posterior density for each of the simulated data 
sets to determine which input model parameters are consistent with 
the observations. We can construct frequentist confidence intervals 
based on Monte Carlo simulations (i.e., comparing the posterior 
distributions for the synthetic data sets to each other). 

The basic approach of the T analysis is similar to likelihood- 
free methods more commonly used in approximate Bayesian com- 
putation. In this case, we do have a likelihood which allows us to 
sample from the posterior probability distribution using standard 
MCMC. We compare the posterior densities calculated for several 
simulated data sets to the posterior density of the actual observa- 
tions, so as to assess the accuracy and bias of the standard MCMC 
analysis. 

The problem of biased eccentricity estimators for nearly cir- 
cular o rbits is familiar from prev ious studies of binary stars. In par- 
ticular. lLucv & Sweeney! d 197 ll) investigated the possibility of mis- 
takenly assigning an eccentric orbit to a circular spectroscopic bi- 
nary due to inevitable measurement errors. As many spectroscopic 
binaries may have been affected by tidal circularization, they sug- 
gested assigning a circular orbit to any system for which the eccen- 
tricity credible interval contained 0. When studying a population of 
systems for which circular orbits are common, this approach sig- 
nificantly reduces the chance of erroneously concluding the system 
has a non-zero eccentricity. One obvious disadvantage of this ap- 
proach is that it would lead to a negative bias for systems where the 
eccentricity is of order cr/K, where a is the typical measurement 
precision and K is the velocity amplitude. For binary stars, cr/K 
may be small enough that this is not a significant concern. For exo- 
planets, where cr/K may be as small as ~ 2-3, such a procedure 
would result in a significant negative bias for many systems. The T 
analysis offers an alternative approach, which may be particularly 
useful when analyzing the eccentricity distribution of a population 
of planetary systems. 

For th e sake of comparison, w e consider a modernized ver- 
sion of the ILucv & Sweenevl dl97 ll ) approach which is based on 
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Figure 5. Cumulative distributions functions (CDFs) of eccentricities from 
different methods. The solid red line is for ejjcMC. adopted from MCMC 
method (Appendix IAU . The dotted red line is similar to MCMC, but any 
eccentricity with e\ mver of for a 95% credible interval is assigned to 0. 
The blue lines are for T analysis, where the solid line is for er from the T 
analysis (Appendix IAU and the dotted line is similar, but any eccentricity 
with eimter of for a 95% confidence interval is assigned to 0. 

the posterior distribution from a standard MCMC analysis or the 
confidence interval from our T analysis. We construct a histogram 
or cumulative distribution of the eccentricities for a population of 
systems, using a single summary statistic for each system: the pos- 
terior mean for the standard MCMC a nalysis or the eccentricity that 
maximizes the T statistic. Following ILucv & Sweenevl ( 1 197 lh . we 
adopt an eccentricity of zero for any system for which the 95% sig- 
nificance level (r method) or the 95% credible interval (MCMC) 
includes e = 0. The cumulative distribution functions of the eccen- 
tricities usi ng different methods are plotted in Fig. [5] Based on the 
generalized ILucv & Sweenevl d 197 lh approach. ~81% (70%) of the 
short-period planet systems in our sample are consistent with circu- 
lar o rbits using the T analysis (standard MCMC analysis). Clearly, 
the ILucv & Sweenevl d 197 ll) approach results in a large fraction 
being assigned a circular orbit, largely due to the choice of a 95% 
threshold. The fraction assigned a circular orbit is sensitive to the 
size of the credible interval used when deciding whether to set each 
eccentricity to zero. There is no strong justification for the choice of 
the 95% threshold (as opposed to 68% or 99.9% threshold) and tun- 
ing the threshold to agree with other meth ods negates the primary 
advantage of the ILucv & Sweenevl dl97 lh method, that it requires 
no addit ional computation s . Ther efore, we do not recommend us- 
ing the ILucv & Sweenevl dl97 ll) approach to learn about the ec- 
centricity distribution for a population when cr/K is not large. 



4 TIDAL INTERACTION BETWEEN STAR AND 
PLANET 

Several factors affect the eccentricity distribution of short-period 
planets including tidal interaction between host star and planet and 
possible perturbation of an undetected companion. We will discuss 
how these factors affect the eccentricity distribution and whether 
the effect is observable. 

In order to understand the influence of tidal interaction on ec- 
centricity distribution, we first divide our sample into two subsets, 
one subset contains systems with a long-term RV trend while the 
other subset contains systems that do not show a long-term RV 
trend (Fig. [6). The systems that are noted with a long-term ve- 
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Figure 6. Distribution of short-period single-planet systems in (e,T age lT c i rc ) 
space. Filled circles are systems showing no linear RV trend and open cir- 
cles are systems showing long-term linear RV trends. Different colors indi- 
cate different combinations of Q' t and Q' 



locity trend include 51 PEG, BD -10 3166, GJ 436, GJ 86, HD 
107148, HD 1 18203, HD 149143, HD 68988, HD 7924, HD 99492 
and r Boo. We further divide the no-trend subset into two groups, 
one group is distinguished by T age lT circ > 1, and the other group 
is distinguished by T c , ge /T cjrc < 1, where T cjrc is tidal circulariza- 
tion time scale and T age is the age of the host star. We investigate 
whether there is a significant difference in the eccentricity distri- 
bution between these two groups as expected if tidal interaction 
is an important facto r in sh aping eccentricity distribution. Follow- 
ing Ma tsumura et al U2008h . we estimate Tare using: 



2 Q' M„ ( a 



81 n M,\R P 



Q'\M tJ , 



F. + F„ 



(5) 



where the subscripts p and * denote planet and star, M is mass, R 
is radius, a is semi-major axis, Q' is modified ti dal quality factor 
and n - [G(M, + M p )/a i ] 1/2 is the mean motion. iMatsumura et al.l 
d2008l) adopted 10 6 as a typical value for Q' t for short period plan- 
etary system host stars and considered Q ranging from 10 s to 10 9 . 
We use a = [10 6 , 10 7 ] and Q' p = [10 s , 10 7 , 10 9 ] in our analysis. 
The factors F, and F p are defined in the following two equations: 
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where CI is rotational frequency. For short-period planets one could 
set Qpjot/n = 1 based on the assumption that all the pl anets in our 
samp le have been synchronized since t !Y „ c/ , ~ 10~ 3 t c ,, t iRasio et"al] 
1 19961) . In order to check whether our conclusion is sensitive to 
the choice of £2» iTO ,/n, we conduct calculations with other fl, trot /n 
values in which we choose stellar rotation period to be 3, 30, 
and 70 days for all the stars. We find that this range for n» ro ,/M 
does not change the conclusions in the paper. Therefore, for fu- 
ture discussion, we adopt Q« iTO ,/n = 0.67, which results in stel- 
lar ro tation periods consiste nt with typical values from 3 to 70 
days jMatsumura et alj|2008l) . The uncertainties in £2« iTO ,/n are ac- 
counted for by our subsequent data analysis (Equation. I lOt. And f\ 
and fi are approximated by the equations: 



f l( e 2 ) = (l + ^e 2 + ^eU^)Kl-e 2 )«> 2 , 



Planet radius R p is estimated based on iFortnev et ail d2007h . 
We assume that the planet and host star are formed at the same 
epoch. We assume that planet structure is similar to Jupiter with 
a core mass fraction of 25M e /M,, = 7.86%. Radii of GJ 436 b 
and HD 149026 b are adopted from reference papers b ecause there 
is a f actor of >2 difference between observed values dTorres et al .1 
120081) and theoretically calculated values. Stellar radius and age 
estimations are obtained from the fo llowing sources with descend- 
ing priority: 1, iTakeda et al] ( 120071) : 2, nsted.ipac.caltech.edu; 3, 
exoplanet.eu. The calculated T circ values are presented in Appendix 
IA2l in addition to the results of MCMC analysis of individual sys- 
tem and other stellar and planetary properties. 

We use the eccentricity posterior samples for each system for 
which the standard MCMC analysis was successful (i.e., results of 
§2.1) to construct the eccentricity samples of two groups separated 
by T age /T cjrc . We note that there are considerable uncertainties in 
the estimation of T age and T„ ro so r agf ,/r arc . > 1 does not necessar- 
ily mean that the actual system age is larger than the actual circu- 
larization time. We consider the sensitivity of our results to these 
uncertainties by adopting a probability function: 



P(— ) = 

1 rirr 



1 - 0.5 exp 
0.5 exp 



-T(^-D 



if 321 < 1 



(10) 
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where r\ is a parameter tuning the confidence of T age and T circ esti- 
mation. For example, if T ag e/Tcirc = 2 and rj = 1, then p(T„ ge /T circ ) = 
0.816, meaning there is 81.6% chance that the system is from the 
group of T age /T cjrc > 1 because of the uncertainties in T age and T cjrc 
estimation. Therefore, we take 81.6% of the eccentricity posterior 
samples of the system to construct eccentricity sample for the group 
of T a geI T circ > 1 an d the remaining 18.4% eccentricity posterior 
samples to construct eccentricity sample for group of T age /T C j rc < 1. 
The rj parameter reflects our confidence in r age and T cjrc estimation. 
If we are not very confident in the estimation of T age and T arc , then 
we set 77 to a small value approaching zero, so half of the eccen- 
tricity posterior samples for each system are assigned to the group 
with T age lT c ir C > 1 and the the other half are assigned to the group 
with T age lT circ < 1. After constructing the eccentricity sample for 
the two groups, we use two-sample K-S test to test the null hy- 
pothesis that these two samples from two groups were drawn from 
the same parent distribution. The results (Table [T) show that we 
are unable to exclude the null hypnosis at a low p value (statis- 
tic of two-sample K-S test) because of the small effective sample 
size (7V'~8). If A„ MA = 0.2, where A,,,,,, is the maximum difference 
between cumulative distribution functions of two groups, we can 
exclude the null hypnosis at p = 0.05 only if N' is more than 44. In 
comparison, our current sample size is inadequate to draw a statis- 
tically significant conclusion on whether or not the the groups are 
from the same parent distribution. However, we do see a hint of a 
difference between cumulative distribution functions of two groups 
(Fig. [7] left), there are more systems with low-eccentricity for the 
group with Tage/Tdrc < 1, which is a consequence of tidal circular- 
ization. We also find that the conclusion is unchanged for a wide 
range of tj, Q' t , and Q values. 

We conduct similar test for two subsets distinguished by 
whether or not a long-term velocity trend is recognized and find the 
similar result that our current sample size is inadequate to draw a 
statistically significant conclusion on whether or not the the groups 
are from the same parent distribution. Again, we see a hint of a 
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Figure 7. The left panel compares the cumulative distribution function of 
eccentricity for two groups of planets: 1) Tagelicirc > 1 (solid line) and 2) 
TagelTcirc < 1 (dotted line). The right panel compares cumulative distribu- 
tion functions of eccentricity for two subsets of planets: planets without a 
long-term RV slope (dashed line) and planets with a long-term RV slope 
(dash-dotted line). 



Table 1. Two-sample K-S test result. A„, a> is the maximum difference be- 
tween cumulative distribution functions of eccentricity for two groups sep- 
arated by Tage/Tarc = 1. N' is the effective sample size, calculated by 
(A^i • N2)I(N\ + A^), p is the significance level at which two-sample K-S 
test rejects the null hypothesis that the two eccentricity samples are from 
the same parent distribution; r\ is a parameter tuning the confidence of T age 
and Tcjrc estimation. 
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difference between cumulative distribution functions of two sub- 
sets (Fig. UJ right) although it is not statistically significant. There 
are more systems with low-eccentricity for subsets showing no sign 
of external perturbation. The maximum difference between the cu- 
mulative distribution functions A max is 0.123, N' is 8.37 and K-S 
statistic is 0.999. In that case, we need an effective sample size of 
1 19 in order to make a statistically significant conclusion (p=0.05). 
In other cases, larger sample size is required since A,„ flv is less. 

We have shown that any difference in eccentricity distribu- 
tion depending on expected time scale for tidal circularization or 
the presence of additional bodies capable of exciting inner planet's 
eccentricity is not statistically significant, although this may be a 
consequence of small effective sample size. The data are also con- 



sistent with the argument that both factors play roles in affecting 
the eccentricity distribution. 



5 ECCENTRICITY DISTRIBUTION 

We seek an analytical function that is able to approximate the ob- 
served eccentricity distribution for short period single planetary 
systems in the framework of Bayesian inference. For this purpose, 
we first exclude systems showing long-term RV trends to reduce 
the effect of perturbation on the estimated eccentricity distribution. 
We also assume that the distribution of T age /T cjrc in our sample is 
representative of short-period single-planet systems. 

Using the posterior samples of eccentricity from standard 
MCMC analysis, we obtain an observed eccentricity probability 
density function (pdf) f(e) by summing the posterior distributions 
together. While not statistically rigorous, this provides a simple 
summary of our results. Logarithmic binning is adopted because 
the shape of /(e) at low eccentricity is of particular interest. The 
uncertainty <x(e) for each bin is set by assuming a Poisson distri- 
bution. Then, we use a brute-force Bayesian analysis to find the 
most probable values of parameters for the candidate eccentric- 
ity pdf /'(e) that approximates the observed eccentricity pdf. In 
the observed eccentricity pdf, there is a pile-up in small eccen- 
tricities near zero and a scatter of nonzero eccentricity less than 
0.8. Therefore, we use a mixture of two distributions: an expo- 
nential pdf f ex po(e,A) = (1//1) • exp(-e//i) to represent the pile- 
up of small eccentricity near z ero and either a uniform distribu- 
tion or a Rayleigh distribution jjuric & Tremainell2008r) to repre- 
sent the population with sizable eccentricities. We assume a uni- 
form distribution for parameters in prior 1.(9), where 8 is vec- 
tor containing the parameters for fie). Our results are not sen- 
sitive to the choice of priors. We adopt Poisson likelihood for 
each bin in the form of fp isson( n \ v) = (V ■ exp(-y))/rc!, i.e., 
L(e t \§) = /p o ,i.„„,(/'(e,|0) • AT;/(e,) • AM, where N is total num- 
ber of posterior samples. The value of /'(e,|#) ■ N is rounded if 
it is not an integer. The posterior distribution of 8 is calculated as 



M | „ Id 



et\6)Y*$ , where Mis 



the number of bins. We have explored a range of bin size from 
10 to 40 bins in logarithmal space and conclude that the results of 
Bayesian analysis do not change significantly with choice of bin 
size. Since the plausible values of parameters for f'(e) are limited, 
we do not have to explore a large parameter space. Therefore, a 
brute-force Bayesian analysis is practical. 

We apply Bayesian analysis to three different planet popu- 
lations for different 77 values (^=0.001,1,1000) assuming Q„ = 
10 7 and Q' p = 10 7 : 1) systems without a long-term RV slope 
and Tage/Tdrc > 1; 2) systems without long-term RV slope and 
tageharc < U 3) the union of 1) and 2). The results of Bayesian 
analysis are presented in Table [2] Fig [8] shows marginalized proba- 
bility density of different parameters for analytical eccentricity dis- 
tribution with a mixture of exponential and Rayleigh distributions 
(/'(c) - a ' fexpa( e >X) + (y ~ a )'fmyi( e > cr ))- a is the fraction of expo- 
nential distribution component and f ray i(e, cr) represents Rayleigh 
distribution with the form of f ray i(e,cr) = (el a 2 ) ■ exp (-e 2 /2cr 2 ). 
Group 1 (potentially circularized system) is represented by dashed 
lines of different colors indicating different values of 77 while group 
2 (systems unlikely to have been circularized) is represented by col- 
ored solid lines. Group 1 and group 2 are well mixed if t/=0.001, 
i.e., a loose constraint is applied on the boundary of T age /T circ = 1. 
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Figure 8. Marginalized probability density functions of parameters for ana- 
lytical eccentricity distribution with a mixture of exponential and Rayleigh 
pdfs. Uncircularized systems (group 2) are represented by colored solid 
lines while circularized systems (group 1) are represented by colored 
dashed lines. Different colors indicate different 77 values, red: 77 = 0.001; 
green: 77 = 1; blue: 77 = 1000. Solid black lines are for union of group 1 and 
2. 



Figure 9. Marginalized probability density functions of parameters for an- 
alytical eccentricity distribution with a mixture of exponential and uniform 
pdf. Uncircularized systems (group 2) are represented by colored solid lines 
while circularized systems (group 1) are represented by colored dashed 
lines, red: 7/ = 0.001; green: 77 = 1; blue: 77 = 1000. Different colors in- 
dicate different 77 values. Solid black lines are for union of group 1 and 2. 



When adopting 77=0.001, the marginalized posterior pdfs for both 
groups 1 and 2 approach the marginalized pdf for group 3 (solid 
black lines; the union of group 1 and 2, including all planets with- 
out a velocity slope). When we adopt larger 77 values (77= 1 or 1000), 
a, the fraction of the exponential component of the pdf is greater for 
group 1 (Fig[8]top: blue and green dashed lines) than group 2 (Fig 
[8] top: blue and green solid lines). This is consistent with the hy- 
pothesis that significant tidal circularization affected group 1. For 
group 1, the fraction of the exponential component of the pdf is 
consistent with unity, implying that the eccentricity distribution for 
planets from group 1 can be described by an exponential pdf. In 
comparison, there is a substantial fraction of Rayleigh component 
(40%) for the analytical function describing eccentricity distribu- 
tion for uncircularized systems (i.e., group 2). Similar conclusions 
are also drawn for the analytical eccentricity distribution with a 
mixture of exponential and uniform distribution (Fig[9) in the form 
of /'(«) = a ■ f expo (e,A) + (1 - a) ■ f mif (e,/3), where /„ m /(e,j8) is 
uniform distribution with lower boundary of and upper boundary 
of p. Fig. [TO] shows the cumulative distribution of the eccentric- 
ities sample based on summing the posterior eccentricity samples 
of each system (green solid line) and the CDFs of the two analytic 
functions (blue dotted and dashed) with the parameters that maxi- 
mize the posterior probability (see Table|2] subset 3). For compar- 
ison, the cumulative distribution of e MCM c from standard MCMC 
analysis (Appendix IAU is shown in red. The difference between 
the red and the green line cannot be distinguished at 0.05 signifi- 
cant level in a K-S test (N'=46). 

In order to check whether the analytical function with the most 
probable parameters is an acceptable approximation to the eccen- 
tricity distribution for short-period single-planet systems, we gen- 
erate test samples from the analytical functions and compare the 
resulting eccentricity samples to observations. For each test sample 
there are N' eccentricities following the distribution of the analyti- 
cal function f'(e\9), where N' is the effective sample size. Each ec- 
centricity is perturbed by an simulated measurement error that fol- 
lows the distribution of posterior samples for our analysis of the ac- 
tual observations. The test sample is then compared to the observed 
eccentricity samples obtained by standard MCMC analysis using 
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Figure 10. Cumulative distributions functions (cdf) of eccentricities from 
different methods. MCMC: most probable eccentricities, &mcmc > adopted 
from MCMC method (Appendix |AT}; MCMC, sum: eccentricities by sum- 
ming up posterior distribution samples of each system; AF1 : cdf of the ana- 
lytical function with the most probable parameters, f'(e) = a ■ f eX po(fi, X) + 
(\-a)- f ray i(e, cr); AF2: cdf of the analytical function with the most probable 
parameters, f(e) = a ■ f„ po (e, A) + (1 - a) ■ f lmi f(e,P). 

two-sample K-S test. We report in Table [2] the percentage of tri- 
als in which eccentricity samples generated by our analytical func- 
tion can not be differentiated from the observed eccentricity sample 
at 0.05 confidence level. All the candidate analytical functions we 
have tested are able to reproduce the observed eccentricity distribu- 
tion in more than 99.6% of the trials at 0.05 confidence level. We 
conclude that the best-fit analytical function is an adequate approx- 
imation to observed eccentricity distributio n. We also compare to 
analytical eccentricity distribution used in IShen & Turner! j2008l) 
/'(e) ~ [1/(1 + e)" - e/2 a ] in which a=4, although it is not specifi- 
cally for short period single planetary systems. Similar to what we 
did in previous test, we found that in 39.0% of the tests, the eccen- 
tricity samples generated by analytical functions can not be differ- 
entiated from the observed eccentricity sample at 0.05 confidence 
level. 
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Table 2. Bayesian analysis results. Group 1: systems without long-term RV slope and T age /T c j rc > 1; group 2: systems without long-term RV slope and 
Tage/Tcirc < U group 3: union of 1 and 2. Numbers in bracket are uncertainties of the last two digits. The fraction column reports the percentage of trials in 
which eccentricity testing samples generated by analytical function can not be differentiated from the observed eccentricity sample at 0.05 confidence level. 
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Figure 11. Distribution of short-period single-planet systems in our sample 
in period-eccentricity space (top: p < 10 day; bottom: p > 10 day). Systems 
that are not consistent with zero eccentricity according to T analysis (filled 
circles) or MCMC analysis (open circles) are marked with corresponding 
names. We also include transiting planets (marked as cross) for comparison. 
Transiting systems not consistent with zero eccentricity are marked with 
numbers: 1, WASP-18; 2, WASP-12; 3, WASP- 14; 4, HAT-P-13; 5, WASP- 
10; 6, XO-3;7,WASP-6;8,WASP-17;9, CoRoT-5; 10,HAT-P-11; ll.HAT- 
P-2. 



From the results of Bayesian analysis, there is a clear differ- 
ence in a, the fraction of exponential distribution, between group 1 
and 2, suggesting the role played by tidal circularization. Group 1 
with T aS elTcirc > 1 shows more planets with near-zero eccentricities 
(a ~ 75%) as compared to group 2 (a ~ 55%) with T„ ge /T circ < 1 
(Table |2j. Since the eccentricity samples tested were perturbed by 
measurement errors, the analytical function we found can be in- 
terpreted as an approximation of the underlying eccentricity distri- 
bution of short-period single-planet systems. However, the actual 
parameter values and uncertainties in the analytical function are de- 
pendent upon the quality of observation and the number of systems 
in the sample of short-period single planets. As the measurement 
precision and the sample size improve, we will be able to better 
constrain the values of parameters in the analytical function which 
approximates the underlying eccentricity distribution. 



6 DISCUSSION 

The median eccentricity of short-period single-planet systems in 
our sample is 0.088. When compared to median eccentricity for all 
the detected exoplanets 0.15, it suggests that the population may 
be affected by tidal circularization. We use eccentricity estimated 



by I" analysis in the following discussion since it is a less biased 
method for accessing small eccentricity. Fig. QT] shows the period- 
eccentricity correlation. We would expect planets with sufficient 
short period to be tidally circularized. While this is generally true, 
there are 3 (17.6%) planets with P < 4d and non-circular orbits: 
GJ 436, HD 149026 and r Boo. For GJ 436, it is sus pected that an 
outer companion may b e pumping the eccentricity dManess et al.l 
l2007l : lRibasetalJl2008l) . Similarly, observations of t Boo b are in- 
consistent with circular orbit, but might be explained by the pertur- 
bation ofanjmserarompanion indicated by a long-term RV linear 
trend (Butler et al. 2006J). Secondary eclipse timing indicates that 
the e ccentricity of HD 149 026 is quite small but inconsistent with 
zero dKnutsonetaT1l2009l) . Considering planets with orbital peri- 
ods up to 10 days, there are 4 additional systems that are not cir- 
cularized. HD 118203, HD 68988 and HD 1852 69 might be due to 
perturbations by additional bodies in the system l lButler et al .120061) 
while the non-zero eccentricity of HD 162020 b may be attributed 
to a different formation mechanism l lUdrv et al.ll2002l) . With period 
longer than 10 days, there are 8 (44.4%) planets with T age /T circ > 1 
which have non-zero eccentricity and no detected long-term linear 
trend. In comparison, there are 8 (28.6%) eccentric planets with or- 
bital periods less than 10 days. The increasing fraction of recogniz- 
ably eccentric orbits as period increases is suggestive of decreasing 
tidal circularization effect, but large sample of planets is required to 
draw firm conclusion. The discovery of over 705 planets candidates 
by the Kepler mission presents an excellent opportun ity to ana- 
lyze the eccentricity distrib ution of short-period planets jFord et all 
l2008t iBorucki et alj|201 j) . We briefly consider transiting planets. 
We note that 11 of 58 (19.0%) transiting systems as of June 2010 
are not consistent with zero eccentricity, and the orbital periods for 
all transiting planets but 3 (i.e. CoRoT-9b, HAT-P-13c HD 80606b) 
are less than 10 days. We infer that the tidal circularization pro- 
cess might be effective for isolated planets with orbital period of 
less than 10 days. Alternatively, the planet formation and migra- 
tion processes for short-period giant planet may naturally lead to a 
significant fraction of nearly circular orbits, even before tidal effect 
takes place. 



It is worth noting that HD 17156 b dFischer et aT]|2007l) has 
one of the most eccentric orbits among short-period planet sys- 
tems in spite t he fact that T aee /T circ could be well over 10 (Fig [6] 
red). However, iBarbieri et al. I l2003) found no indications of addi- 
tional companions based on o bservations of directing imagi ng, RV 
and astrometry measurement. lAnglada-Escude et al .1 ( 120101) inves- 
tigated the possibility that a 2:1 resonant orbit can be hidden by 
an eccentric orbital solution. It is interesting to explore such pos- 
sibility on this particular system to solve the discrepancy of high 
eccentricity and r age to T circ ratio. Another possibility is that the 
system is in the process of circularization that began well after the 
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star and planet formed (e.g., due to planets scattering). However, 
we are cautious in drawing conclusions since T age lT circ could be 
less than unity (Fig[6] green and blue). 



7 CONCLUSION 

We apply standard MCMC analysis for 50 short-period single- 
planet systems and construct a catalog of orbital solutions for these 
planetary systems (Appendix I A2t . We find general agreement be- 
tween MCMC analysis and previous references with the primary 
exception being cases where eccentricity was held fixed in previ- 
ous analysis. We develop a new method to test the significance of 
non-circular orbits (r analysis), which is better suited to perform- 
ing population analysis. We find the eccentricity estimations from 
T analysis are consistent with results from both standard MCMC 
analysis and previous references. 

Our results suggest that both tidal interactions and external 
perturbations may play roles in shaping the eccentricity distribu- 
tion of short-period single-planet systems but large sample sizes 
are needed to provide sufficient sensitivity to make these trends sta- 
tistically significant. We identify two analytical functions that ap- 
proximate the underlying eccentricity distribution: 1) mixture of an 
exponential distribution and a uniform distribution and 2) a mixture 
of an exponential distribution and a Rayleigh distribution. We use 
Bayesian analysis to find the most probable values of parameters 
for the analytical functions given the observed eccentricities (Table 
O. The analytical functions can be interpreted as the underlying 
distribution of eccentricities for short-period single-planet systems. 
Thus, the analytical functions can be used in the future theoretical 
works or as priors for eccentricity distribution. 
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